Imaging the water content of the inner 


astronomical units of HL Tau 


Stefano Facchini", Leonardo Testi?, Elizabeth Humphreys?^^, 


Mathieu Vander Donckt?, Andrea Isella’, Ramon Wrzosek”, 
Alain Baudry, Malcom D. Gray?, Anita M. S. Richards!?, 
Wouter Vlemmmings! 


“Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 
16, Milano, Italy. 
?Dipartimento di Fisica e Astronomia “Augusto Righi", Università di 
Bologna, Viale Berti Pichat 6/2, Bologna, Italy. 
? European Southern Observatory, Karl-Schwarzschild-Strasse 2, 
Garching bei München, D-85748, Germany. 
^ Joint ALMA Observatory, Alonso de Cordova 3107, Santiago, Chile. 
? European Southern Observatory (ESO) Vitacura, Alonso de Cordova 
3107, Santiago, Chile. 

6 Space sciences, Technologies & Astrophysics Research (STAR) 
Institute,University of Liége, Liége, Belgium. 
"Department of Physics and Astronomy, Rice University, 6100 Main 
Street, MS-108, Houston, 77005, TX, USA. 

*Laboratoire d'Astrophysique de Bordeaux, Univ. de Bordeaux, CNRS, 
B18N, allée Geoffroy Saint-Hilaire, Pessac, 33615, France. 

? National Astronomical Research Institute of Thailand, 260 Moo 4, T. 
Donkaew, A. Maerim, Chiangmai, 50180, Thailand. 
19JBCA, University of Manchester, M13 9PL, UK. 

H Department of Space, Earth and Environment, Chalmers University of 
'Technology, Góteborg, SE-412 96, Sweden. 


*Corresponding author(s). E-mail(s): stefano.facchini@unimi it; 


28 


29 


30 


31 


32 


33 


34 


35 


36 


37 


38 


39 


40 


41 


42 


43 


44 


45 


46 


47 


48 


49 


50 


51 


52 


53 


54 


55 


56 


57 


58 


Abstract 


The water molecule is a key ingredient in the formation of planetary systems, 
with the water snowline being a favorable location for the growth of massive 
planetary cores. New ALMA data of the ringed protoplanetary disk orbiting the 
young star HL Tauri show centrally peaked, bright emission from water vapour 
in three distinct transitions of the main water isotopologue for the first time in a 
quiescent planet forming disk. The spatially and spectrally resolved water content 
probes gas in a thermal range down to the water sublimation temperature. Our 
analysis implies a stringent lower limit of 3.7 Earth oceans of water vapour 
available within the inner 17 astronomical units. We show that, due to the high 
dust column density and absorption, our observations probe the water content in 
the atmosphere of the disk, indicating the main water isotopologue as the best 
tracer to unveil water vapour in planet forming regions. 


The water molecule is undoubtedly one of the most important molecular species in 
the whole universe. Being an extremely efficient solvent, water had a key role in the 
emergence of life as we know it on our planet. For this reason, the chemical charac- 
terization of exoplanetary atmospheres is often focused on detecting this particular 
molecular species [1-3]. Formed by the common H and O atoms, water plays a funda- 
mental role in the physics of the formation of planetary systems, due to its very high 
abundance in both gaseous and icy forms [4, 5]. Theoretical models predict that at the 
location of the phase transition from gaseous to solid form, dust grains can accumulate 
and grow very efficiently, promoting the fast formation of planetary cores. Across this 
particular radial location, called ‘snowline’, grains can drastically change their drift 
and fragmentation velocity, composition, and opacity. In synergy with vapour radial 
diffusion [6], these physical discontinuities can lead to the accumulation and growth 
of dust grains into planetesimals [7, 8]. The position of the snowline also defines the 
chemistry of the available planet building blocks. Since the H20 molecule is the major 
elemental oxygen carrier in the disk, its desorption and freezing affect the elemental 
C/O ratio in both the gas and solid phases [9-11]. 

Because of its large binding energy, the H2O transition from ice to gas happens a 


few astronomical units (au) from the young star where the midplane temperatures are 
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in the range from 100 to 200 K, making it the last major ice component to sublimate. 
However, the proximity to the host star makes the detection of the snowline compli- 
cated even in the closest star forming regions. Both cold and warm water lines have 
been detected in a few disks by Herschel (see [12, 13] and references therein), Spitzer 
[14], JWST [15-17] and ground based observatories [18], but the low angular resolu- 
tion did not allow robust inferences about the extent of the water snowline. Observing 
directly water emission from the ground is complicated by the high water vapour 
content of the Earth’s atmosphere, resulting in strong telluric absorption. To circum- 
vent this problem, most programs at mm wavelengths have focused on attempting the 
detection of the rarer H}8O and HDO isotopologues, leading to the clear detection 
of spatially resolved water isotopologue emission in the outbursting source V883 Ori 
[19]. In quiescent sources, except the candidate detections in the AS 205N and HL 
Tau disks, no detection of thermal emission at (sub-mm) wavelengths has been 
reported in the literature [20-23]. 

In this paper we focus on the text-book case of HL Tau, the first protoplanetary 
disk imaged at very high angular resolution (~ 0.025") with the Atacama Large Mil- 
limeter/submillimeter Array (ALMA) [24]. The disk shows a spectacular pattern of 


concentric rings. With the source being young (S 1 Myr), and the dynamical stel- 


lar mass being relatively high (2.1 + 0.2 Mo, [25]), the inner disk temperatures are 
expected to be warm, due to irradiation and accretion heating. Warm and hot water 
has been detected in HL Tau both by Herschel in the Far Infrared (FIR, [26]) and 
by ground-based high resolution spectroscopy in the Mid-Infrared (MIR, [27]), with 
the lines not being spatially resolved. With this article, we report the detection 
of three rotational water lines in the inner regions of the HL Tau protoplanetary disk 
obtained with ALMA. The lines are spectrally resolved. Analysis of the interferomet- 


ric data confirms that the extent of the water emission is confined within a prominent 
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gap seen in the HL Tau continuum intensity. These new data present the first spa- 
tially resolved images of the emission from the main water isotopologue (H20) in a 
protoplanetary disk, and pave the way to a new observational strategy to characterize 


the water vapour content of terrestrial planet forming regions. 


Observations 


We observed HL Tau in two different ALMA bands (Band 5, originally developed with 
the goal of studying water in the local Universe [28], and Band 7), to target three 
transitions of water (two lines of para-water, and one line of ortho-water): p-H2O 
313 — 220 and 515 — 42», at 183.31 and 325.15 GHz, respectively, and o-H3O 1059 — 936 
at 321.22 GHz. The first two lines are expected to trace warm water vapour 
outside the water snowline, while the third line is predicted to detect hot 
water within the water snowline [29, 30]. We also observed a rotational transition 
of the water isotopologue p-H}80 at 322.46 GHz. The molecular coefficients of the 
lines, and sensitivity of the observations, are reported in Table 1 (and Table $2). 
After self-calibration, we imaged both the continuum and the continuum- 
subtracted lines with the CASA software [31]. The continuum images are shown in 
Fig. 1. All the H20 lines and the H480 line were imaged with CASA 6.2.1 with natural 
weighting, to maximise point source sensitivity. The 183 GHz line presents a synthe- 
sized beam of 0.500" x 0.442" (PA —3.0deg), and was imaged with a channel width 
of 0.8 km s^. The resulting beam for the 325 GHz water line is 0.640" x 0.491" (PA 
—42.1 deg), with a 1 km s^! channel spacing. For the high excitation 321 GHz line, only 
one long baselines execution block was available, and the beam is 0.067" x 0.061" (PA 
12.2 deg), with a 5 kms"! channel spacing. The H}80 line was imaged with several 
different channel spacings. To have a one-to-one comparison with the main isotopo- 
logue line, in the paper we show the results with a channel width of 1 kms-!. The 


resulting beam is 0.779" x 0.626" (PA —46.7 deg). 
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'The spectrum of the two lowest excitation lines was extracted over a circular area 
with 0.7" radius and is shown in Fig. 2. Both the 183 GHz and the 325 GHz lines 
are clearly detected across multiple channels, with the lines being centered on the 
systemic velocity of the system (7.1kms~!; [32, 33]). The 325 GHz line shows an 


l. as seen for other lines at similar upper energy 


absorption component at ~ 10kms~ 
levels [33]. The 183 GHz line shows a width of ~ 12 kms", while the higher frequency 
lines show a slightly broader emission. The 183 GHz line spectrum exhibits a peak 
signal-to-noise ratio (snr) of ~ 9.6, with an rms of 13.2 mJy over the flux density 
extraction area. The 325 GHz line spectrum instead shows a peak snr of — 5.8, with 
an rms of 14.0 mJy. The higher energy line at 321 GHz does not show a detection when 
integrating over the same area, due to the much smaller beam and the resulting high 
rms. We thus extracted a spectrum over a circular area with radius 0.06". The peak 
snr is ~ 4.1, with an rms of 3.0 mJy. This transition shows the broadest line profile 
among the three detected lines, consistent with originating from the inner regions with 
higher Keplerian velocities. In all cases, the integrated intensity map shows a strong 
detection with a peak co-located with the dust continuum peak (see Figure 1). We 
obtain a peak snr of 21.4, 19.8 and 8.1 in the integrated intensity maps of the 
183, 325 and 321 GHz lines, respectively. The line fluxes extracted over a circular area 
centered on the continuum peak are reported in Table 1. 

The intensity weighted velocity maps (moment one maps) of the 183 and 325 GHz 
lines are shown in Fig. 1. In the high snr map of the 183 GHz line, disk rotation is clearly 
detected, with position angle and systemic velocity in agreement with other brighter 
lines from the same dataset [33], indicating that a displacement of the photocenter of 
blue-shifted and red-shifted channels is detected at the current resolution. 

For the H1*O line in Band 7, we extracted the spectrum with the same methodology 
described for the three main isotopologue lines, over a 0.7" radius circular area. A 
moment 0 map was computed over the same spectral range as the main isotopologue 


line. The line was not detected. 
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The spatial distribution of water vapour 


The integrated intensity morphology of the 183 and 321 GHz lines are remarkably dif- 
ferent (see Figure 1), but so are the angular resolutions of the two moment maps. In 
order to derive radial profiles of the respective integrated intensities, we used two dif- 
ferent approaches. First, we focused on the highest snr line of the sample (at 183 GHz), 
which is also the coldest and therefore expected to trace the largest spatial extent. We 
averaged in frequency the interferometric data of the continuum-subtracted line in the 
same frequency range used to compute the moment 0 map. We then fitted the line inte- 
grated intensity visibilities with a simple Gaussian model using the galario package 
[34], after fixing the inclination and position angle to 46.7 deg and 138.0 deg, respec- 


tively, as derived from high angular resolution continuum observations [24]. The fit 


converges well to a Gaussian with og = 0.12 0.01" (see Fig. 3; og is the standard 


deviation of the Gaussian function). At a distance of 140 pc, this corresponds to 


oq = 16.8 + 1.4au. For the 321 GHz line, the much higher angular resolution allowed 
us to compute the radial profile of the integrated intensity by azimuthally averag- 
ing the moment zero map, after de-projecting it by the known inclination (using the 
GoFish package [35]). 

Figure 3 compares the two integrated intensity profiles with the Band 7 
azimuthally averaged continuum intensity profile. Both line profiles are clearly cen- 
trally peaked, indicating that the water vapour emission must originate above the 
optically thick continuum from the disk midplane. The lower excitation line is signif- 
icantly more extended, showing that the water emission has a radially decreasing 
temperature (excitation), and that the high excitation line is not optically thick out- 
side the central beam. The same line shows detectable emission out to 0.3" when 
boosting the snr by azimuthally averaging. The 183 GHz temperature gradient is con- 
firmed by fitting the high snr spectrum with a Keplerian disk model with a brightness 


temperature gradient, assuming that the line is close to being optically thick. 


166 


167 


168 


169 


170 


171 


172 


173 


174 


175 


176 


177 


178 


179 


180 


181 


182 


183 


184 


185 


186 


187 


188 


189 


190 


191 


Declining temperature profiles are preferred to flat ones, in agreement with the derived 
integrated intensity profile. 

The warm brightness temperature profile of the 183 GHz line, which far exceeds 
the midplane temperatures obtained by analysis of multiwavelength continuum data 
[36], indicates that the water vapour we are tracing originates in the warm disk upper 
layers. This is further supported by the peak of the 321 GHz emission, which is slightly 
shifted from the continuum peak (see Figure 1, bottom right panel). Even though 
the two are consistent within the astrometric precision of the data, the apparent shift 
agrees with tracing water vapour on the side of the disk closer to the observer, which 
is in the North-East [25], unocculted by the optically thick dust continuum. These 
findings set a stringent upper limit on the radius of the water snowline at 
17au. A more accurate determination will require simultaneous forward- 
modelling of the radial and spectral profiles of all three (sub-)mm lines 


with the aid of thermo-chemical codes. 


Water column density 


From the three main water lines, we computed the rotational diagram, under the 
assumption of optically thin emission and uniform excitation temperature across the 
energy range. While the very inner regions of the line emission are likely opaque, the 
bulk of the emission from the three lines cannot be optically thick, since this would 
imply a flux ratio that is equal to the square of the frequency ratio (in Rayleigh-Jeans 
regime and Local Thermodynamical Equilibrium — LTE), which we do not observe. 
While masing cannot be excluded to partly contribute to the emission, in the high 
spectral resolution spectra of the 183 and 325 GHz lines we have not identified indi- 
vidual narrow spectral features which are in general a good signature for maser action 
together with high flux density. From the high densities of the inner disk, masing is not 


expected for the three lines analyzed here, and the only contribution could originate in 
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the very upper layers at low volume densities where collisional quenching of the 
masing action is less probable. In the rotational diagram, we used degeneracy 
quantum numbers that account for a 3:1 ortho-to-para ratio, and a parti- 
tion function that considers all states within the same assumption [37]. No 
rescaling of the ortho- and para- line fluxes is needed to compute the total 
water column. We accounted for 1096 absolute flux systematic uncertainties in the 
line fluxes. 

The rotational diagram does not provide a unique solution for the rotation tem- 
perature, as shown in Figure 4. The Monte Carlo Markov Chain (MCMC) exploration 
individuates two distinct temperatures, which realistically indicate a continuous gra- 
dient in the excitation temperature of the water vapour. Fitting either the two lower 
energy lines or the two higher energy lines separately, the rotational diagram indi- 
cates excitation temperatures of 914712 and 7901127 K, respectively, in line with the 
two classes of solutions obtained in the joint fit. The lower temperature solution is 
sensitive to colder water vapour in the range of expected desorption temperatures of 
water ice in space, suggesting that the bulk of the water emission from the 183 GHz 
line traces water gas in the proximity of the snow surface. The higher temperature 
solution is driven by warm gas in the upper layers of the terrestrial planet forming 
regions of the disk, which are well imaged by the high resolution 321 GHz line inte- 
grated intensity map, and likely by the lines being close to be optically thick. The 
excitation temperature of the warm gas is in broad agreement with MIR line 
fluxes in the 12.37 — 12.41 um range on the same source (in particular the 
o-H5O 16413 — 15114 line with Eup = 4948 K [27]). These high temperatures 
allow for water vapour formed in situ via gas-phase reactions [38, 39]. 

All the rotational diagrams robustly constrain the column density of water gas (in 
the optically thin limit and above the optically thick continuuum). The joint 


fit shows a total water column density logio (Nain /cm?) = 16.414696 within 0.7" 
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(~ 100au) from the star. This corresponds to ~ 3.7 Earth oceans (7.1 x 10? 
lunar masses) of water vapour. Given the optically thin assumption, this has to be 
considered as a tight lower limit. Since most of the emission originates from < 0.12" 
(17 au), assuming that the entirety of it is confined within this radial range we obtain 
an averaged column density of logio (Nui /cm?) ~ 18.10*005, when accounting for 
the disk inclination. 

The non-detection of the H15O line can determine an upper limit of the aver- 
age optical depth expected for the HzO 325 GHz line (7y,0). Assuming that TH O = 
530 Ty180 (using the oxygen isotope ratio measured in the solar wind [40]), we obtain 
that Tmo < 14. By turning the argument around, the optically thin assumption for 
the 325 GHz line provides an upper limit of N(H30)/N(H1*O)- 40 in HL Tau, well 
in agreement with oxygen fractionation levels in our own solar system [41]. The non 
detection of H15O shows that observational campaigns aimed at targeting water in 
inner disks of quiescent disks with ALMA should privilege the main isotopologue as 
a first choice, with follow-up observations of robust detections targetting 
HDO and H}80 to derive more accurate column densities and set stringent 
limits on the water deuteration [19]. 

Assuming that the bulk of the water emission originates from < 17au from the 
star, we can compare the total mass of water (My,0) to the mass of dust Maus 
estimated from multi-wavelength continuum analysis of ALMA and VLA data [36]. 
By using the dust surface density from this study, we obtain Maus ~ 13MEartn, which 
leads to a water-to-dust mass ratio of My,0/Maust ~ 6 x 107°. This number is much 
lower than the expected water abundance in inner disks (water-to-dust mass ratio 
~ 10-?). The optical thickness of water can only marginally alleviate the problem, 
given the non-detection of H4°O. Continuum subtraction could also marginally reduce 
the water brightness temperature [42], but in this case is not expected to reduce the 


water fluxes by more than a factor of 2 (see Figure 3). The low water-to-dust ratio 
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further strengthens the interpretation that with ALMA we are probing only the upper 
layers of the disk, above the optically thick screen of the dust continuum emission, 
which shows optical depths between ~ 5 — 10 at 0.9 mm within the inner 17 au [36]. 
The large cavity seen in the HDO and H3O integrated intensity maps of 
V883 Ori [19] further supports that the observation of a large column of 


water is hindered by optically thick continuum in the inner disk. 


Conclusions 


These new ALMA data reveal high significance detections of three distinct rotational 
transitions of the main isotopologue of water vapour in the inner regions of the ringed 
HL Tau disk. These observations pave the way to the characterization of the water 
content of the inner regions of protoplanetary disks. The tremendous angular resolu- 
tion and sensitivity of the ALMA telescope, even in spectral ranges of low atmospheric 
transmission, are providing the first spatially and spectrally resolved images of the 
vapour of the main water isotopologue in a planet forming disk. Analysis of the mor- 
phology of the water emission, of the spectrum of the highest snr line, and of the 
excitation conditions, jointly indicate that the (sub-)mm lines are probing warm gas in 
the disk upper layers above the water snow surface, with a radially decreasing temper- 
ature profile. The non detection of H15O and the low water-to-dust ratio in the inner 
17 au show that the observations are probing marginally optically thick gas above the 
opaque dust continuum emission. These results highlight that the water content of 
quiescent protoplanetary disks at (sub-)mm wavelengths is most efficiently unveiled 
by targeting the main isotopologue, in particular in disks with high continuum optical 


depths within the water snowline. 
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Transition Vo Eu Ch. width rms^ Flux Mask rad? 


(GHz) (K) (kns ^!) (mJybm !) (mJykms '!) (") 
p-H2O 313-229 183.31009 204.7 0.8 6.90 973 + 89 0.70 
p-H2O 515-499 325.15290 469.9 1.0 3.71 1332 + 89° 0.70 
0-H20 1029-936  321.22569 1861.3 5.0 0.96 679 + 135 0.28 
p-HIFO 5,5-45, 322.46517 467.9 10 1.75 < 33.67 0.70 


Table 1 Observed H20 isotopologue transitions. The upper state energies are taken from [43]. The 
notation for the water energy levels in the vibrational ground state is Jz, jy. Notes: ^ Obtained 
over one single channel. ? Circular radius used to extract the line flux. ^ This flux measurement 
corrects for the absorption identified at ~ 10 kms-!. The flux obtained without accounting for 
absorption is 929 + 89 mJy km s^ !. 4 3c upper limit. 
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Fig. 1 Top. Left: 1.7mm continuum image of HL Tau. Center: integrated intensity map of the 
183 GHz water line. Right: intensity weighted velocity map of the 183 GHz water line, after a 40 
clipping on individual channels, where disk rotation is clearly detected. Center. Same as top panels, 
for the 0.94mm continuum and the 325 GHz water line. The intensity weighted velocity map in 
this case is computed after a 30 clipping. Bottom. Left and center: same as top panels, for the 
0.94 mm continuum and the 321 GHz water line. No moment one map is shown, due to low snr. Right: 
zoom-in of continuum intensity, with [4,5,6,7,8]o contours of the 321 GHz line moment 0 map, with 
c = 13.3mJy beam! km s^ !. The rms associated to the the integrated intensity maps of 


the 183 and 325 GHz lines are respectively: 28.2 and 46.3 mJy beam"! km s^ !. 
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Fig. 2 Top left: spectrum of the 183 GHz water line, extracted over a circle with radius of 0.7” 
centered on the continuum peak. Top right: spectrum of the 325 GHz water line, extracted over the 
same area, highlighting the mirrored (flipped) version of the spectrum with the dashed-dotted line. 
Bottom: spectrum of the 321 GHz water line extracted over a circle with radius of 0.06" centered on 
the continuum peak. The velocity range on the r—axis is different in the bottom panel. The grey 
dashed line in all panels shows the systemic velocity of 7.1 kms-!. In the top left of each panel 
2c scale bars are reported for reference. 
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Fig. 3 Left: integrated brightness temperature (Tj) radial profile of the 321 GHz line, recon- 
structed integrated Tj radial profile of the 183 GHz line, and Tj, profile of the 0.94 mm continuum 
emission. The thick orange line shows the best fit model assuming a Gaussian integrated intensity 
profile. Thin lines show randomly sampled realizations of the posterior distribution. The lines in 
the top right show the beam major axis. For the 183 GHz line it portrays the smallest spa- 
tial scales to which the uv-plane analysis is sensitive, which are — 2.3 smaller than the 
major axis of respective natural beam. Right: Re-centered and de-projected visibilities of the 
integrated intensity of the 183 GHz water line. The thick red line shows the best fit model reproducing 


the profile shown on the left panel. 
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Fig. 4 Left: rotational diagram of the three water lines, with line fluxes extracted as in the main text. 
The fit does not lead to a unique solution, indicating that the assumption of uniform temperature is 
inadequate. Right: posterior distribution of the rotational diagram fit . The dashed lines indicate the 
16th, 50th and 84th percentiles of the marginalized posterior distributions. 
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Methods 


Observations, data reduction and imaging 


HL Tau was observed in both Band 5 and Band 7 with the ALMA Program 
2017.1.01178.S (PI: Humphreys), targeting the two para-water lines at 183.31004 GHz 
and 325.15297 GHz, respectively. HL Tau was also observed in band 7 to target the lat- 
ter water line with Program 2022.1.00905.S (PI: Facchini), together with the H480 line 
at 322.46517 GHz (see Table S1). The molecular coefficients for the three transitions 
are reported in Table 1 and $2. 

Within the 2017.1.01178.5 ALMA program, the source was observed in Band 5 on 
September 21, 2018, for a total integration time of 46 min, with 43 antennas and base- 
lines ranging between 15 m and 1.4 km. The weather conditions during the observation 
were exceptional, with a median precipitable water vapour (PWV) column during 
the observations of ~ 0.19 mm. J0423-0120 was used as amplitude and bandpass cal- 
ibrator, whereas J0510--1800 for phase referencing. The Band 5 spectral setup had 6 
spectral windows (spws), five of which targeting different molecular transitions, includ- 
ing the 183 GHz water line, and a 1.875 GHz-wide spw for continuum observations at 
170.004 GHz. Within the same program, HL Tau was observed in Band 7 in two spec- 
tral setups. The first one on August 12, 2019, with a time on-source of 31 min, with 48 
antennas and baselines ranging between 41 m and 3.6 km. With a median PWV column 
of 0.4 mm, J0431+1731 was used to cross-calibrate phases, and J0538-4405 for ampli- 
tude and spectral response. This first spectral setup consisted of four spws in FDM 
mode; three of them with a maximum bandwidth of 1.875 GHz, and one with 1920 
244 kHz channels, targeting the water 325 GHz line. The second spectral setup was 
used within the same program with an execution block observed on November 24th, 
2017. The median PWV was 0.5 mm. J0519-4546 was used as amplitude and bandpass 
calibrator, and J0440+1437 as phase calibrator. The observation spent 33 min on the 


science target with 49 antennas, with a maximum baseline of 8500 m. The spectral 
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setup consisted of 4 spws in FDM mode, three of them with maximum bandwidth, 
and one of them with 1920 244 kHz channels, targeting the water 321 GHz line. 

Finally, new data were taken in October 2022 with a more compact array in Band 
7 with Program 2022.1.00905.S, with baselines ranging between 15 and 500 m, and 
a total on-source time of 100 min with 41/45 antennas (in two execution blocks). 
The median PWV column was ~ 0.3. J0423-0120 was used for flux and bandpass 
calibration, while J0431+1731 was used for phase referencing. With four spws, one 
of them was centered on the 325.15 GHz water line, while another spw targetted the 
H150 line at 322.46517 GHz. 

The Cycle 6 (9) data were calibrated by the ALMA pipeline using CASA v5.4 
(v6.4) [31]. For the band 7 data, we combined the data from the two cycles. For both 
bands, we first self-calibrated each of the three execution blocks in phase, combining 
all spectral windows after flagging spectral ranges associated to line emission, and 
combining all scans. Following [44], we then aligned the data by fitting a Gaussian 
to the continuum, and shifting the phase-center to the continuum maximum with the 
fixvis and fixplanet tasks. We then self-calibrated the two short-baselines (in the 
case of band 5, only one set of baselines is available) execution blocks in both phase 
and amplitude, reaching a peak snr of 9950 (a — 50096 improvement). In the case of 
band 7, we then combined the long-baseline execution block from the Cycle 6 data, 
and self-calibrated the data again in both phase and amplitude, reaching a peak snr of 
5300 (a 21096 improvement). We took particular care with the amplitude calibration, 
where the gain solution with scan-length intervals greatly improved the data quality. 
The models for self-calibration were constructed with CLEAN with Briggs weighting 
(robust=0.5). The gain solutions were then applied to the full spectral data. 

The Band 5 continuum data were imaged with robust—0.0. With a synthesized 
beam of 0.364" x 0.312" (PA —9.5deg), the Band 5 (1.70mm) continuum presents 


an rms of 33 uJy, with a peak snr of 2457. The band 7 (0.94 mm) data were imaged 
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with robust=—1.0. The data exhibit an rms of 36 uJy, with a peak snr of 403, over 
a synthesized beam of 0.037" x 0.029" (PA 1.5 deg). A flux density for both images 
was obtained over an elliptical area with a semimajor axis of 1.3", and semiminor axis 
and position angle (PA) to trace the disk inclination and PA (46.7 and 138.0 deg, 


respectively; [24]). Without accounting for absolute flux calibration uncertainties, we 


obtain a flux density of 323.0 + 1.4 mJy and 1677.7 + 0.6 mJy at 1.70 and 0.94 mm, 
respectively, where the uncertainties have been computed as standard deviations of 
randomly selected masks with the same area of flux density extraction over emission- 
free regions of the continuum maps. 

While the water lines have been imaged with natural weighting (see main text), 
several additional attempts with a range of uv-tapers were performed to increase the 
sensitivity to extended emission, but they did not show any feature undetected with 
natural weighting. Both integrated intensity (moment zero) and intensity weighted 
velocity (moment one) maps were generated for the water lines. Moment zero maps 


! without any clip- 


were computed by integrating channels between —2 and 16 kms^ 
ping (Fig. 1) for the Band 5 line, between —6.4 and 20.6kms"! for the 325 GHz 
line, and between —10.4 and 24.6 kms^! for the 321 GHz line. For the two brighter 


lines, we integrated the moment zero map over a circular area centered over the emis- 


sion peak and a radius of 0.7". We obtain line fluxes of 973 + 89 mJy kms"! and 
929 + 89 mJy km s^! (for the 183 and 325 GHz lines, respectively). Since the 325 GHz 


line shows an absorption red-shifted component, we also computed the underlying 


flux by considering the blue-shifted side only, and multiplying it by a factor of two. 


The resulting flux is 1332 + 89 mJy kms^!. The uncertainties on the line fluxes were 
computed by bootstrapping over 100 circular apertures in emission free regions of the 
map, and do not account for absolute flux calibration uncertainties. The same oper- 


ation was applied to the 321 GHz line, using a smaller 0.28" radius extraction area. 
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The resulting flux is 679 + 135 mJy kms^!. The flux is much lower than the tenta- 
tive detection by [20] with the Submillimeter Array (SMA), where however the line is 
shifted by 30km s^! from the systemic velocity, indicating that the proposed emission 
may be originating from a large scale flow which we filter out in our high resolution 
data. For the 183 and 325 GHz lines, intensity weighted velocity maps were generated 
using the bettermoments package [45] after applying 4 and 3c clipping to individual 
channels, respectively. 


The non detection of the H50 322 GHz line is shown in Figure 5. 


Fitting of the 183 GHz line spectrum 


Given the high snr of the 183 GHz line, we fitted its spectrum by analytically comput- 
ing model spectra of a geometrically thin Keplerian disk, similarly to [46]. To do 
So, we assumed that the peak brightness temperature of the line decreases with radius 


as a power-law: 


T(R) = To (=) B (1) 


For a given radius R, and azimuth $, we assumed that the line intensity follows a 


Gaussian distribution velocity: 


B,(T) umg(v — vk proj)” 
I = l 2 
(R, 4,0) = F erp (8 e), (2) 
where B,(T) is the Planck function at temperature T, kp is the Boltzmann 
constant, d is the distance of HL Tau (140 pc), ymy is the mass of the water molecule, 


and UK, proj can be written as follows: 


GMo\ UA 
we = (SRE) ico (3) 
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where i is the source inclination (46.7 deg). We considered an optically thick 
limit when assuming a thermal broadening with a kinetic temperature 
equalling the brightness temperature. The flux density of the line can then be 


computed at every velocity v by integrating across the whole disk: 


27 Rout 
F(v) = cosi f f I(R,,v)RdRdọ, (4) 
0 Rin 


where cosi accounts for the geometrical projection on the sky. We fixed Rin 
to 0.1 au (but the model is not sensitive to this value for reasonably small radii), and 
sampled the disk with 150 points in radius, and 550 in azimuth. We then convolved 
the models with a Gaussian kernel with the same channel width as in the data, and 
sampled them at the same velocities. We kept three free parameters in the fitting 
procedure: To, q and Rout. We fitted the spectrum shown in Fig. 2 with the emcee 
package, using flat priors on the three free parameters: [10,1500] K, [0,3], [2,200] au, 
respectively. We used 30 walkers, 1000 steps of burn-in, and 1000 additional steps to 
sample the posterior distribution. Fig. 6 shows the best fit model, and 100 random 
draws extracted from the posterior distribution. While the fit does not manage to 
constrain the outer radius of the emission, we obtain To = parse and q = 0.92+9:30. 
The fit highlights a negative brightness temperature gradient in the radial profile, as 
seen in the reconstructed integrated intensity profile (see Figure 3), and as hinted by 


the rotational diagram shown in Fig. 4. 


Parametric fit of the integrated intensity profile in the visibility 


plane 


In order to extract the radial extent of the 183 GHz line, we performed a paramet- 
ric fit of its integrated intensity radial profile in the visibility plane, by exploiting 


the galario package. After averaging the continuum-subtracted visibilities in the 
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same spectral range used to compute the moment zero map, we fitted the visibility 
data by Fourier transforming a projected integrated intensity profile in the same uv- 
points sampled during the observations. We modelled the radial profile with a simple 
Gaussian prescription: 

R? 

J(R) = Jo exp (-zx) , (5) 
where we considered four free parameters: Jo, oc, and the disk center (ARA and 
ADec). We fixed the inclination and position angle to the ones obtained from high 
resolution continuum imaging [24]. The fit was performed with the emcee package, 
where the Jo parameter was sampled in log-space. We used the following flat priors 
on the parameters: log;9(Jo/steradian) € [1,20], og € [0,1.5]", ARA € ([—0.4, 0.4]", 
ADec € [-0.4, 0.4]". The posterior distribution was sampled with 50 walkers and 1000 
steps, after 1000 steps of burn-in. The MCMC exploration of the posterior space well 


converges, as shown in Fig. 7. 


Rotational diagram and optical depth constraints 


'To compute the rotational diagram of the water molecule, we use the same approach as 
in e.g., [47, 48]. In the optically thin assumption, we can compute the column density 


Nihin and the rotational temperature T ot by measuring the integrated flux S, Av: 


Nihin = 


4m S, Av Q(Trot) Ey 
Ayhe Q Gu EP Trot ' 


(6) 
where Q is the solid angle used for flux extraction (see previous section), Ay is the 
Einstein coefficient of the considered transition, h and c are the Planck constant and 
the speed of light in vacuum, Q is the partition function, Eu is the upper state energy 


(in K) and gu the upper state degeneracy. Using the relation between Nuthin and 


Ninin, the same equation can be written in logarithmic form [49]: 


Nuthin 
In 


g — ln Nihin 2 In Q(Tiox) = Eu/Trot- (7) 
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We performed a linear regression using the emcee sampler [50] in the 
[In (Nu tnin/gu), Eu] space to extract Minin and Trot: We used the molecular coeffi- 
cients reported in Table S3. In the fitter, the partition function was determined 
with a cubic spline interpolation across the rotational temperatures listed in [37]. 
The same approach was used for the rotation diagram analysis with six to 
eight transitions of o- and p-water in evolved stars [51]. In the Monte Carlo 
Markov Chain (MCMC) sampling, we used 128 walkers, and 2000 steps (after 1000 
steps of burn-in). While Figure 4 shows the result and the marginalized posterior dis- 
tributions of the fit of all three lines, Figure 8 portrays the individual fits on the two 
colder and warmer lines, respectively. 

To compute the upper limit on the optical depth of the 325 GHz water line, we 
exploited the non-detection of the H15O line, which has almost identical molecular 
coefficients. Assuming that the column density of the main isopotologue line (74,0) 


is equal to 530 x rgiso [40], we can write: 


(S, ^v)m,o s 1—e 77H20 
(SA)upo ^ 1— e740 780" 


(8) 
Using the 30 upper limit on the H480 transition, and the measured flux of the H2O 
line (after correcting for absorption, since the absorbing column of H1*5O will have the 
same scaling factor of 530), a numerical solution of the equations leads to TH,o < 14, 


as reported in the main text. 
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» Extended data 


——1 km s™ channel 
—— 10 km s^! channel 
==- Vsys 


Flux Density (Jy) 


-20 -10 0 10 20 30 
Velocity (km/s) 


Fig. 5 Spectrum of the 322 GHz H480 line, extracted over a circle with radius of 0.7” centered on 
the continuum peak, as for Fig. 2. Two spectra are shown, from channel maps with two different 
channel widths. The grey dashed line indicates the systemic velocity of 7.1 km s^! [32, 33]. In the 
top left 2c scale bars are reported for reference. 
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Fig. 6 Left: spectrum of the 183 GHz H20 line, as in Fig. 2. The uncertainty associated to each data 
point is shown in the left part of the spectrum as an errorbar. A random sampling of 100 profiles from 
the posterior distribution of the MCMC fit is shown, with the dark red line indicating the bestfit 
model. Right: marginalized posterior distribution of the fitted parameters, where no constraint can 
be retrieved on Rout from this analysis. 
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Fig. 7 Marginalized posterior distribution of the galario fit of the visibility points of the integrated 
intensity of the 183 GHz line (see Fig. 3). 
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Rotational diagram for two colder lines 
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Rotational diagram for two warmer lines 
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Fig. 8 Rotational diagrams and marginalized posterior distribution of the MCMC exploration for 
individual fits on the two lowest and two highest energy lines, respectively. 
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Program ID Lines Int. time PWV . Bp/flux cal. Phase cal. Max. bl 


(min) (mm) (m) 
2017.1.01178.8 — p-H30 313—220 46 0.2 J0423-0120 . J05104-1800 1398 
2017.1.01178.8 — p-H30 515-425 31 0.4 J0538-4405 . J04314-1731 3637 
2017.1.01178.8 0-H2O 1059-936 33 0.5 J0519-4546 . J04404-1437 8547 
2022.1.00905.8 p-H20 515—422 100 0.3 J0423-0120 J0431+1731 500 


p-HP50 515—422 


Table 2 Observations IDs and Execution Blocks properties. 


Transition du Au Q(100K) Q(200K) Q(300K) Q(400K) 
(s!) (K) 
p-H20 315-2529 7  3.59x 10-6 35.1 97.4 178.1 274.6 


p-H30 515-422 11 1.15 x 1075 E - E _ 
o-H3O 1029-936 63 6.21 x 10-9 , - 7 g 
p-H}8O 515-459. 11 1.05 x 10-5 - i - - 


Table 3 Molecular coefficients of the observed H20 isotopologue transitions are from the JPL 
database [43], with radiative coefficients from [52] and the updated partition function 
Q from the ExoMol database [37], which well agrees with the one by [52] in the 
temperature range explored in this paper. Only some representative values are reported in 
this table. The degeneracy quantum number of the ortho-state assumes an 
ortho-to-para ratio of 3. The partition function is computed over all states with the 
same ortho-to-para ratio. 


Declarations 


Data availability 

Al the ALMA data are publicly available on the ALMA archive 
(https://almascience.nrao.edu/aq/). 

Code availability 


The python packages used in the data analysis are all publicly available. The 


calibration and fitting scripts can be obtained by SF upon reasonable requests. 
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